#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>

#define RANDZEROONE drand48()
#define RANDONEONE ((2.0f*drand48())-1.0f)
#define RANDZEROFIVE (drand48()-0.5f)

#define V(side) ((side)*(side)*(side))

enum RAND { ZEROONE, ONEONE, ZEROFIVE };

//-------------------------------------------------------------------------------------------------------
// Variable Host Declaration
//-------------------------------------------------------------------------------------------------------
int init_type_j;

int nspin;
int side;
int nsteps;
int overrelaxiter;
int outfreq;
int dumpfreq;
int seed;

int ns_j;

size_t ms_s;
size_t ms_j;

REAL* Ra;
REAL* Rb;
REAL* Rc;
REAL* Ba;
REAL* Bb;
REAL* Bc;

REAL* Jpx;
REAL* Jpy;
REAL* Jpz;
REAL* Jmx;
REAL* Jmy;
REAL* Jmz;

REAL* aRand; 
REAL* bRand; 
REAL* cRand; 
REAL* Rand;

int sdbL;
int sdbLs;
int sdbLm1;
int sdbLsm1;
int maskL;
int maskL2;
int maskL3;
int Ldb2;
int Lsdb2;

REAL *Energy;
REALENE  totene;

REAL h;
REAL Beta;
REAL zero, one, two, three, zerofive, zerotwentyfive;
REALENE zeroEne; //MESSA SOLO PER IL CALCOLO DELL'ENERGIA IN DOUBLE COME PROVA

bool onlyOverrelaxation = FALSE;
bool onlyHeatBath = FALSE;

int randNum; // DA TOGLIERE MESSO PER TEST RIMETTERE LOCALE

int verbose=0;



//-------------------------------------------------------------------------------------------------------
// Function Declaration
//-------------------------------------------------------------------------------------------------------
void OverrelaxRed(REAL *a, REAL *b, REAL *c, REAL *an, REAL *bn, REAL *cn, REAL *Jpx, REAL *Jpy, REAL *Jpz, REAL *Jmx, REAL *Jmy, REAL *Jmz);
void OverrelaxBlue(REAL *a, REAL *b, REAL *c, REAL *an, REAL *bn, REAL *cn, REAL *Jpx, REAL *Jpy, REAL *Jpz, REAL *Jmx, REAL *Jmy, REAL *Jmz);

void HeatbathRed(REAL *a, REAL *b, REAL *c, REAL *an, REAL *bn, REAL *cn, REAL *Jpx, REAL *Jpy, REAL *Jpz, REAL *Jmx, REAL *Jmy, REAL *Jmz);
void HeatbathBlue(REAL *a, REAL *b, REAL *c, REAL *an, REAL *bn, REAL *cn, REAL *Jpx, REAL *Jpy, REAL *Jpz, REAL *Jmx, REAL *Jmy, REAL *Jmz);

void loopXRed(REAL *a, REAL *b, REAL *c,
              REAL *an, REAL *bn, REAL *cn,
              REAL *Jpx, REAL *Jpy, REAL *Jpz,
              REAL *Jmx, REAL *Jmy, REAL *Jmz,
              REAL *t_a, REAL *t_b, REAL *t_c,
              int istart, int iend, int yzoff,
              int spxoff, int smxoff, int spyoff, int smyoff, int spzoff, int smzoff);

void loopXBlue(REAL *a, REAL *b, REAL *c,
	       REAL *an, REAL *bn, REAL *cn,
	       REAL *Jpx, REAL *Jpy, REAL *Jpz,
	       REAL *Jmx, REAL *Jmy, REAL *Jmz,
	       REAL *t_a, REAL *t_b, REAL *t_c,
	       int istart, int iend, int yzoff,
	       int spxoff, int smxoff, int spyoff, int smyoff, int spzoff, int smzoff);

void loopXRedHB(REAL *a, REAL *b, REAL *c,
		REAL *an, REAL *bn, REAL *cn,
		REAL *Jpx, REAL *Jpy, REAL *Jpz,
		REAL *Jmx, REAL *Jmy, REAL *Jmz,
		REAL *t_a, REAL *t_b, REAL *t_c,
		int istart, int iend, int yzoff,
		int spxoff, int smxoff, int spyoff, int smyoff, int spzoff, int smzoff);

void loopXBlueHB(REAL *a, REAL *b, REAL *c,
		 REAL *an, REAL *bn, REAL *cn,
		 REAL *Jpx, REAL *Jpy, REAL *Jpz,
		 REAL *Jmx, REAL *Jmy, REAL *Jmz,
		 REAL *t_a, REAL *t_b, REAL *t_c,
		 int istart, int iend, int yzoff,
		 int spxoff, int smxoff, int spyoff, int smyoff, int spzoff, int smzoff);

void ComputeHEnergy(REAL *a, REAL *b, REAL *c, REAL *an, REAL *bn, REAL *cn,
		    REAL *Jpx, REAL *Jpy, REAL *Jpz,
		    REAL *Jmx, REAL *Jmy, REAL *Jmz,
		    REAL *energy,
		    int even);


void Init(void);
int findlog2(int n);
void FreeHost(void);
void Initialize_Spin(REAL* sa_a, REAL* sa_b, REAL* sa_c);
void Initialize_J(void);
REAL rand_gauss(void);
void Parameters_Reading(char *inputfile);
void Calculate_Sizes(void);
void Dump_Spin(int count);
void ComputeRandom(REAL* rand_ptr, int mode);


